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Abstract. We present a non-Markovian quantum trajectory method for treating 
Ft , , atoms radiating into a reservoir with a non-flat density of states. The results of an 



o\ 



example numerical simulation of the case where the free space modes of the reservoir 
are altered by the presence of a cavity are presented and compared with those of an 
extended system approach. 
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c-| ■ 1. Introduction 

The radiative properties of an atom depend critically on the nature of the 
electromagnetic field surrounding the atom. It is therefore of interest to consider atoms 



situated in cavities |I| 0, [3|, |], |5], wave guides and optical fibers || [7] or photonic 
band gap materials [§, P, [lO], 11, |l2| where the atoms have been shown to experience 
a very different process of decay from that in free space. In many of these cases the 

rapidly varying mode structure of the radiation reservoir invalidates the Born-Markov 

c3 ' 



approximations usually employed to treat atoms radiating into free space [[0J. As a 
consequence of this problem, many theoretical works have considered these situations 



by treating the complete unitary dynamics of atom(s) and field [T(| |TT|, [12]. 

When the Born-Markov approximation applies, the reservoir degrees of freedom 
can be eliminated to derive a master equation for the reduced dynamics of the atom 



T4|l . Quantum trajectories |L6|] (stochastic Schrodinger equations [[OJ or Monte Carlo 



wavefunctions [T7^) have been found very successful in numerically solving this master 



equation |L9, 20, 21, 22]. A generalization of the quantum trajectory method beyond 



the Born-Markov regime has been made independently by Diosi et al. E3[] and by 



the present authors [24] . This generalization opens the door to a treatment of atomic 
radiation into general reservoirs in terms of a reduced atomic system. By adhering 
to the measurement interpretation of quantum trajectories it is possible to bypass the 
intermediate step of deriving a master equation and go directly to a non-Markovian 
equation for the state of the system conditioned on a continuous measurement process 
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24| . Where the Born approximation does not hold this state is no longer pure owing to 



the presence of an entanglement between the reservoir and the atom. 

In this paper we summarize the derivation of the non-Markovian trajectories and 
explicitly demonstrate how to simulate the trajectories for an atom radiating into the 
above mentioned structured continuum. The simulation method is viable if the number 
of time steps (where the time step At is small on the time scale of the damping) per 
memory time, N = T m /At, is not too large, as keeping track of the system state a long 
time in the past is computationally demanding. 

2. Quantum trajectories from a measurement perspective 

In the rotating wave approximation the coupling of a two level atom to an 
electromagnetic field can be described by the interaction picture Hamiltonian 

fl>(t) = in[^(t)(r / (t)-e(t)«T}(t)] l (i) 

where o~i(t) is the interaction picture atomic lowering operator and 

m=T,9Mt )e-^-^ (2) 

k 

is referred to as the driving field fl25| . The kth mode of the field fy., where [&&, b' k ,} = Skk'i 



is coupled to the atom with a strength g^ and has a frequency of oscillation u>k- 

Quantum trajectories can be derived from a measurement theory perspective by 
considering continuous measurements made on the electromagnetic field that has been 



emitted irretrievably from the atom p^] . The probability of getting a set of results 
R(t,t Q ) over the interval [to,t] given a set of inputs 7(t,to) is, 

Pri(Mo) = M*o)|«L(Mo)nRi(Mo)Mto)), (3) 

where 

n M (Mo) = (R(Mo)|£/mt(Mo)|l(Mo)> (4) 

is the evolution operator for the atomic system conditioned on the particular results 
and inputs, and U{ n t(t, to) is the unitary evolution operator of the total system of atom 
and reservoir. 

Consider a process in discrete time where the time between any two points on our 
grid is At. We can simulate a continuous measurement process by determining the 
conditional probability of getting a particular result at time t given a set of previous 
results. This conditional probability can be written as 

Pai(t,to) = (rm°(t)} 

P R i(t-At,t ) (^ (t-At)|^ (t-At)>' 
where the unnormalized wave function satisfies the interaction picture equation 

$ (t)) = $°(t - At)) + AQ m (t, t - T m ) o $(t - T m )>, (6) 



Non-Markovian quantum trajectories 3 

where we have assumed that there is a finite time in the past, T m , before which the 
state of the system does not effect the probabilities in the present. We will discuss this 
in more detail below and also illuminate the meaning of the superscript. 

The explicit form of ([]) for a photon counting measurement situation with vacuum 
input (where the time step is chosen so that the probability of two particles being 
detected during one time step is negligible) is 

$°{t)) = $°{t - At)) - Ata / (t) t f ds f m (t - s)T | a I (s)l[A(t J )V 1 o ${t - T m )> 



+AN{t) 



(7) 



J^jts h {t _ S ) T I a I (s)l[A(t ] )V I o $(t - T m )) - \4>°(t - At)) 

where the interaction picture operators inside the curly brackets, T{- • •}, are time- 
ordered. The restriction AiV(t) 2 = AN(t) describes a point process in which at time t 
either one photon is detected, AN(t) = 1, or no photons are detected, AN(t) = 0. The 
tj are the discrete times of previous detections during the time interval [t — T m ,t) and 

A(tj) = J* ds h(tj - s)a I (s), (8) 



t-T m 



where 



h(t-s) = Y,9ke- Wk{t - s \ (9) 

k 

is the single particle Green's function for a photon propagating from the atom to the 

detector. We call this quantity the response as it represents the weighting of a particular 

photon emission time given an instantaneous measurement at time t. The interaction 

picture evolution operator has been reduced to the form 

Vo = exp \ - / dsi / ds 2 oj(si)/m(si - s 2 )a I (s 2 ) \ , (10) 

I Jt-Tjn Jt-T m J 

where 

Ut-s) = Y,\9k\ 2 e- tuJk{t - s \ (11) 

k 

is the single particle Green's function for a photon that is emitted and then reabsorbed 



again by the system. This we call the memory function of the system |25 . 

A trajectory is simulated by calculating the probability of a click of the detector 
(^ (t)|^°(t)) A 7V(t)=i x At (or no click of the detector (^ (t)\^°(t)) AN{t)=0 ), dividing 
by (ip°(t — At)\ifj°(t — At)) and comparing this conditional probability with a random 
number between and 1. If the probability is greater (less) than the random number 
then we say that a photon has been irretrievably emitted from the system at time t. 
This result becomes part of the permanent measurement record and is incorporated 
into the subsequent evolution as an additional factor, A(t). The rate of fall off of the 
response and the memory function determines the time t — T m , where T m is referred to 
as the memory time. We can neglect the effect of earlier emissions on the probabilities 
at time t. The state of the system during the memory time is difficult to define because 
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Non-Markovian quantum trajectories 

subsequent measurements could alter the state. However, the state before time t — T t 
will not change with future measurements so we can fix a state at time t— T m conditioned 
on measurements up until time t. It is only in the Markovian case that the conditioned 
evolution operator Qri factorizes at time t — T m , as in this case the memory function 
decays to zero during a single time step and we can put T m = At. In general, the 
entanglement between the system and bath at the time t—T m means that the conditioned 
state of the reduced system is mixed. The o-product of (f|) represents a summation over 
separate conditioned evolution operators acting on separate kets in the system Hilbert 
space, 

Att m {t,t- T m ) o $(t - T m ) = £ An[ R (t,t - T m )|^(t - T m )>. (12) 

i 

The kets of the system with different superscript labels are orthogonal as the labels 

represent a different state of the field. The conditioned state of the system is a density 

matrix given by 

(t T) E,l^-r m )X^(t-r m )| 

£,$'(* -T m )|#(t-T m )> 

and conditioned expectation values of system quantities are then given by (A(t)) c = 
Tr{p c (t)i}. 

3. Simulating non-Markovian quantum trajectoies 



In the above section (and in our previous work [23]) we have given a formal derivation of 
non-Markovian trajectories. The present task is to determine an algorithm to simulate 
a trajectory numerically. Equation ([/]) clarifies the connection between a measurement 
event and the processes of emission and absorption going on in the atom. At time 
t there is a chance of a photon being detected which corresponds to a photon being 
irretrievably emitted at some time in the past. The times of emission are weighted by 
the response function (see the third term on the RHS of (0)). There is also the chance 
of a photon not being detected which corresponds to a photon being emitted at some 
time in the past and then reabsorbed by the system. The emission times are weighted 
by the memory function (see the second term on the RHS of (0)). 

In a simulation, the measurement and the emission events take place at discrete 
times on a grid and the number of time steps per memory time is given by N = T m /At. 
In figure [l| we have given a diagrammatic representation of the three simplest cases. 
To keep track of whether a photon corresponding to a measurement event has or has 
not been emitted at every discrete time we introduce a notation consisting of a binary 
number associated with each possibility where a 1 denotes that a photon has already 
been emitted and a the opposite. Each place in the binary number corresponds to 
a different measurement event. As an example, we consider the state of the system 
at time t — T m in the N = 3 case. A schematic of the situation is given in figure 
0. There are four possible paths that the system can take at time t — 2 At and a ket 
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is associated with each of these possible paths. This example demonstrates a general 
property of these trajectories, in that, if there are N grid points in a memory time and 
each of these points corresponds to a measurement event then there are 2 Ar_1 possible 
paths for the system to take at the time t — T m . One therefore needs 2 7V_1 kets to 
fully specify the conditioned state. The rapid growth of this quantity with increasing 
iV puts a fundamental limitation on this method as a ket must be stored for each of 
these possibilities. Although the kets are in the Hilbert space of the system the different 
labels correspond to orthogonal states of the reservoir because they either correspond to 
distinct measurement results or to different numbers of photons in the field and therefore 
the complete state is an incoherent sum of the kets and is given by fll3|) . 

The above notation is also useful for determining an algorithm to compute the 
forward propagation of the kets. As an example we consider an explicit calculation of 
the measurement probabilities at time t from the state at t — 3At for the N = 3 case. 
We define kets of the form \ip lmn (t)) where the right superscript label corresponds to 
the measurement event at t — 2 At (with a corresponding emission weighting function 
/i = h) the middle superscript to the measurement event at t — At (with associated 
function f 2 = f m ) and the left superscript to a possible event at time t (with function 
/3 = h or f m ), see figure |l[ A first-order algorithm to calculate numerically the first 
increment of the kets in the Schrodinger picture version of (|7|) is 

|^ 000 (t - 2At)) = (1 - iH At)\4> 00 (t - 3At)), 

|^ 001 (t - 2 At)) = (1 - iH At)\4> 01 (t - 3 At)) + oA(0)|^ 00 (£ - 3At)), 

|^ 011 (t - 2A£)> = (1 - iH At)\^ u {t - 3 At)) (14) 

+a (j 2 (At)\i, 01 (t - 3At)) + 7T(0)|^ 10 (£ - 3A£)> 



where we have considered the averaged functions 

7 3 (t 3 -t)= f dsf^-s) (15) 

Jt-At 

as the rate of change of the memory function may be faster than the time scale of 
interest. We must now rearrange our kets before the next time step as we have reached 
the t — 2 At measurement event. We therefore keep only the kets that are labeled with a 
right superscript of 1 as these represent the paths where an emission corresponding to 
the measurement result at t — 2 At has occurred. The other states become redundant. 
The new states are given by 

|^ 001 (t-2Ai))^ |-0 oo (i-2At)), |^ 011 (t-2At))^ |^ 01 (£-2A£)>, ■■• (16) 
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and $2 — * /i and f$ —* f%- The next increment is then 

|^ 00 (t - At)) = (1 - iH At)\^ 00 (t - 2 At)}, 

|^ 01 (t - At)) = (1 - i# At)|^ 01 (t - 2 At)) + (7iT(0)|^ 00 (t - 2 At)), 

|^ 10 (t - At)) = (1 - i# At)|^ 10 (t - 2At)) + a^(At)|^ 00 (t - 2At)), (17) 

|^ n (t - At)) = (1 - iH At)\4> n {t - 2At)) 

+a (7T(0)|^ 10 (t - 2At)) + J 2 {At)\i, m {t - 2At))) . 

This time the measurement event corresponds to no detection so rearranging the states 
must proceed by 

\ft{t - At)) = |^ 00 (t - At)) + at|^ 01 (t - At)), 

$\t - At)) = |^ 10 (t - At)) + J\ft\t _ At)), 

and J2 — > /i. The final increment is 

|^(t)) = (l-itf At)|^ (t-At)), 

|^(t)) = (1 - iH At)\^(t - At)) + aTi(0M°(t - At)). 

The rearrangement at the end to get a final single ket \tp°(t)) will obviously depend on 
whether f\ corresponds to a detection or no detection. For example, if f\ = f m then 
IV^W) = l^ ,0 (^)) + °" I^CO) an d the conditional probability of no detection at time t 
given the previous measurement record is given by (|5p where \ip°(t — At)) is the ket 
corresponding to the realization of the simulated measurement at time t — At. We have 
written out the explicit procedure for the N = 3 case but the basic structure of the 
general case should be apparent from the above example. An increment in each ket has 
the form 

\${t + At)) = (1 - itf At)|^(t + At)) +aJ27n\r(t)}, (20) 

n 

where the f n are taken from the set of responses and memory functions corresponding 
to the measurement record and are multiplied by the appropriate ket, \ip n (t)). The 
summation is over all the possible paths to \ip l (t + At)). The binary notation can be 
exploited to make this step automatic. The kets are then reduced by referring to the 
measurement record. 

In summary, the system state is propagated forward by following all the possible 
paths that the system can take constrained by the fact that the paths must be consistent 
with the measurement record. The measurement probabilities at time t are determined 
by propagating the system forward from time t — T m conditioned on the previous 
measurement record. The time t — T m is chosen because at this is the latest time 
where the state of the system is independent of the outcome of the measurement at 
time t. The outcome of the (simulated) measurement at time t becomes part of the 
permanent measurement record and the system state at time t — T m can be propagated 
forward At in accordance with this record. 2 N ~ l kets must be stored to fully represent 
the conditioned state of the system at any time. By exploiting a special binary notation 
an algorithm can be determined that takes any memory function and response and the 
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associated memory time and generates a trajectory. The Markovian trajectory is then 
a special case of this algorithm when N — 1. 

4. An atom radiating into a cavity 

The practicality of the method is demonstrated by comparing the results of a non- 
Markovian trajectory with the results of a traditional Markovian trajectory for an 
extended system, using the fact that systems that undergoes non-Markovian decay can 
often be embeded in a larger system that undergo Markovian decay |^6|, f27fl . 

The simple test case we consider here is a driven two level atom coupled to a single 
mode of a one-sided optical cavity, where we assume that the atom only emits into the 
cavity. From a non-Markovian perspective the system is the two level atom given by 
the Hamiltonian (in a frame rotating with the classical driving field) 

#o = A^ + -(cx + cx t ), (21) 

where Au is the detuning from the driving field and Q is the strength of the driving 
field. The memory and response functions are given in terms of Fourier transforms by 

Ut - s) = -3= I ( Kiev,.?. ,.^ - Mt -'\ (22) 



2nJ (f) 2 + (cu-z/) 2 



Ht -s) = J^- / duo- P - e -M*-) (23) 

where 7 is the strength of the coupling to the atom, v is the frequency of the cavity in the 
rotating frame and k is the line-width of the cavity. On the other hand by modeling the 
cavity explicitly one could run a Markovian simulation with the effective Hamiltonian 

H e g = H + va)a + \yf^j{ao^ — a) a) — i—a^a, (24) 

where a is the annihilation operator of the cavity mode and the collapse operator is 
y/^ya. A comparison of the results of these two systems are shown in figure |3| and ||. 
The non-Markovian quantum trajectory determines the state of the atom conditioned 
on measurements of the emitted light. By specifying the state of the atom after all 
the light that the atom has emitted has been measured we gain maximum knowledge 
about the atom. In contrast, in the extended system method the state of the atom, 
calculated by tracing the cavity out of the conditioned state of atom and cavity, implies 
a partial erasing of information about the state of the atom. The ensemble average 
of these two methods is equivalent as the average sums the measurement results over 
all times. In figure [5] and ^ we compare a single quantum trajectory with an extended 
system trajectory. 

5. Summary 

We have presented a generalization of the quantum trajectory method for treating 
an atomic system radiating into a reservoir to the non-Markovian regime. We 
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have summarized the derivation of these trajectories from a continuous measurement 
perspective and explicitly demonstrated the algorithm for numerically simulating a non- 
Markovian quantum trajectory. The results of numerical simulations of a simple test 
example were shown to agree with more traditional methods. The method retains many 
of the familiar concepts of Markovian quantum trajectories and, in fact, contains the 
Markovian method as a special case. In the Markovian case the reservoir enters the 
equations of motion for the reduced system as a single rate which is determined by 
integrating the memory and response functions over all time. The non-Markovian case 
requires a more detailed description of the reservoir and the time dependence of the 
memory and response functions becomes important. 

The method described here promises to have wide applicability to many situations 
where weakly non-Markovian behavior (i.e., not too many time steps during a memory 
time) arises. 
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Figure captions 



Figure 1. Schematic of a non-Markovian trajectory. The grid on the x-axis represents 
the propagation of discrete time from left to right. The curves ending on the axis 
represent photons being emitted and then reabsorbed by the system, whereas, the 
straight lines ending above the axis represent the irretrievable emission of a photon 
from the system. The measurement record fixes the events at the end of the lines but 
photons can be emitted at any time after the beginning of the lines up until the lines 
end. The dotted lines ending at t represent the two possible outcomes of a measurement 
at this time. The situation depicted here is: at time t — At no photon was detected 
and at t — 2At one photon was detected. 



Figure 2. Diagrammatic representation of the four possible combinations of emission 
times that contribute to the two measurement results at t — At and t — 2 At of figure EL 
N = 3. The possible regions of emission are represented by the shading. The darker 
shading represents an overlap of the two emission regions. The possibilities are: 00 - 
both the photon detected at t — 2At and the photon reabsorbed at time t — At have 
not been emitted, 10 - only the photon detected at t— 2At has been emitted, 01 - only 
the photon reabsorbed at t — At has been emitted, 11 - both the photon detected at 
t — 2 At and the photon reabsorbed at time t — At have been emitted. 



Figure 3. The real (upper) and imaginary (lower) parts of the correlation function of 
the system (cr (r)cr) determined by an average over a single trajectory of 10 5 detections 
is plotted in (a) for the parameters £1 = 2j, v = 2j, Au> = and k = IO7. The 
spectrum (the Fourier transform of the correlation function) is plotted in (b). Note 
the asymmetry of the spectrum due to the cavity being detuned v ^ from the central 
peak of the Mollow spectrum. The corresponding quantities from the extended system 
treatment are also given (dashed line). For comparison we have also plotted the Mollow 
spectrum for an atom radiating into free space (dotted line). 



Figure 4. The waiting time distribution for the detections is compared with that 
calculated by the extended system treatment (dashed line). 



Figure 5. A typical run of ten detections of a single non-Markovian quantum 
trajectory. We have plotted the probability of a detection as a function of time for the 
parameters: k = 87, £1 — 2j and Aw = v = 0. The same probability determined from 
the extended system method using the same random numbers is also shown (dashed 
line). 
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Figure 6. Plot of the evolution of the conditioned expectation value of a z over 
the trajectory. The extended system method (ESM) and the non-Markovian quantum 
trajectory (NMQT) calculate different states for the atom. Notice that in the NMQT 
case the atom makes full oscillations and also starts emitting before a detection occurs. 
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Figure. 2. 
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Figure. 3. (b) 
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Figure. 5. 
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